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ABSTRACT 

Detection of the 21-cm signal coming from the epoch of reionization (EoR) is challenging 
especially because, even after removing the foregrounds, the residual Stokes I maps contain 
leakage from polarized emission that can mimic the signal. Here, we discuss the instrumental 
polarization of LOFAR and present realistic simulations of the leakages between Stokes pa¬ 
rameters. From the LOFAR observations of polarized emission in the 3C196 field, we have 
quantified the level of polarization leakage caused by the nominal model beam of LOFAR, 
and compared it with the EoR signal using power spectrum analysis. We found that at 134— 
166 MHz, within the central 4° of the field the (Q, U) —>■ I leakage power is lower than the 
EoR signal at fc < 0.3 Mpc“^. The leakage was found to be localized around a Faraday depth 
of 0, and the rms of the leakage as a fraction of the rms of the polarized emission was shown 
to vary between 0.2-0.3%, both of which could be utilized in the removal of leakage. More¬ 
over, we could define an ‘EoR window’ in terms of the polarization leakage in the cylindrical 
power spectrum above the PSF-induced wedge and below fey ^ 0.5 Mpc“^, and the window 
extended up to fcy ~ 1 Mpc~^ at all k± when 70% of the leakage had been removed. These 
LOFAR results show that even a modest polarimetric calibration over a field of view of < 4° 
in the future arrays like SKA will ensure that the polarization leakage remains well below the 
expected EoR signal at the scales of 0.02-1 Mpc“^. 


1 INTRODUCTION 

Five phases of the large-scale universe are imprinted on Hydro¬ 
gen: (i) the primordial phase before redshift 2 ~ 1100—when the 
universe was a hot, dense plasma—that ended when protons re¬ 
combined with electrons releasing the photons that we detect to¬ 
day as a ~ 2.7 K signal known as the cosmic microwave back¬ 
ground (CMB); (ii) the ‘Dark Ages’ (1100 > 2 > 30) when the 
baryonic universe contained mostly neutral Hydrogen and freely 
moving photons; (iii) the ‘Cosmic Dawn’ (30 ^ 2 > 12) when 


the first structures formed; (iv) the ‘Epoch of Reionization’ (EoR; 
12 ^ 2 > 6.5) when high-energy photons emitted by the first 
sources reionized the Hydrogen in the intergalactic medium; and 
(v) the current phase (2 < 6.5) when almost all Hydrogen in the 
universe are ionized (e.g. |Eurlanetto et al.||2006| |Mellema et al.| 
[20T3l|Zaroubi et al.|2012| l. 

The aforementioned highly uncertain boundaries of the EoR 
have been approximated using indirect probes, e.g. CMB polar¬ 
ization at the high -2 end (e.g. [Page et ar]|2007|) and absorption 
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features in quasar spectra at the low -2 end (e.g. |Fan et al.|2006| ). 
However, the new generation low-frequency, wide-bandwidth ra¬ 
dio interferometers have the potential to directly detect the 21-cm 
radiation emitted by neutral Hydrogen during the EoR, redshifted 
to the wavelengths of around 1.5-3 m (corresponding to 200- 
100 MHz), as a differential brightness with respect to the CMB. 
There are several ongoing and planned experiments to detect the 
EoR signal using radio arrays: Giant Metrewave Radio Telescope 
(GMRtQ Low Frequency Array (L0FAR|^ Murchison Widefield 
Array (MWA|^ Precision Array for Probing the EoR (PAPErQ 
21-cm Array (21CMA|^ and the planned Square Kilometre Array 

(SKAjg 


In order to detect the EoR, the effect of all other signals, e.g. 
the Galactic and extragalactic foregrounds, has to be excluded from 
the observed data; spatial fluctuations of the Galactic foreground 
can be 2-3 orders of magnitude higher than that of the EoR signal 
( [Bemardi et al.|20d^|2010||Pober et al.|2013| l which is around 10 
mK within the redshifts 6-10 at 3' resolution ( jPatil et ^|2014| >. 
However, even after removing the foregrounds with high accuracy 
the system noise after even hundreds of hours of integration will 
be an order of magnitude higher than the signal, thereby forcing us 
to aim for a statistical detection of the signal. One of the methods 
for detecting the EoR signal statistically entails removing the fore¬ 
grounds with high accuracy and then measuring the power spec¬ 
trum of the residual which depends heavily on a proper understand¬ 
ing of the systematic and the random (noise) errors associated with 


the observing instrument and foreground removal (e.g. Dillon et 
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In this paper we address the systematic errors due to polarized 
foregrounds associated with the EoR experiment being conducted 
using LOFAR (the LOFAR-EoR project). After taking out the 
bright extragalactic foreground, i.e. the resolved point sources, the 
Galactic foreground can be removed utilizing the fact that the EoR 
signal has significant correlated stmcture along the frequency— 
or equivalently the redshift—axis while the Galactic diffuse fore¬ 
ground is spectrally smooth in Stokes I. However, the Faraday ro¬ 
tated polarized Galactic foreground is not always smooth along fre¬ 
quency and hence a leakage of the polarized emission into Stokes 
1 might mimic the EoR signal (e.g. jjelic et al.|2010| >. Systematic 
errors can cause this leakage in two different ways: direction in¬ 
dependent (DI) and direction dependent (DD). Non-orthogonal or 
rotated feeds of an antenna of an interferometer can cause Q to 
leak into I and vice versa while cross-talk between two feeds can 
cause mixing between all 4 Stokes parameters. As these are DI 
errors, they can be corrected with high accuracy using traditional 
self-calibration. However, the DD errors (DDE) caused by the time- 
frequency-baseline dependent primary beams cannot be corrected 
so easily. In the latter case, an ellipticity of the beam can cause 
I -H- Q mixing while cross-polarization between two orthogonal 
components of the beam can mix all Stokes parameters. 

jCarozzi & Woan| ( |2009^ calculated a full polarization Mueller 
matrix to account for the look-direction dependent polarization 
aberration inherent in a dipole interferometer due to the fact that 


^ http://ginrt.ncra.tifr.res.in/ 

^ http://www.lofar.org/ 

^ http://www.mwatelescope.org/ 
^ http://eor.berkeley.edu/ 

® http://21cma.bao.ac.cn 
® http://www.skatelescope.org/ 


a source sees different projections of a dipole at different times. 
jJelic et aTjpOlO) used this Mueller matrix to calculate the amount 
of leakage to be expected over the field of view of LOFAR and 
found that the leakage should be 0.1-0.7% at 138 MHz within a 
5° X 5° patch of sky around the zenith and should increase to 2- 
20% for an elevation of 45°. If the polarized intensity is ~ IK, 
then a 1.5% leakage would give a polarized emission of ~ 15 mK 
in Stokes I which is comparable to the EoR signal. [Moore et al.| 
( [2013^ simulated the sky with randomly generated Faraday rotated, 
polarized point sources and found that the power of Q —>■ / leakage 
due to the model beam of PAPER that has a FWHM of around 45° 
at 150 MHz is of the order of thousands of [mK]^ which is several 
orders of magnitude higher than the expected EoR signal power. 
Their result turned out to be pessimistic because of their choice of 
the model; in reality, point sources are much more weakly polarized 
at low frequencies ( [Bernardi et al.|20T3) >. 

Here, we predict the level of polarization leakage to be ex¬ 
pected in the 3C196 window of the LOEAR-EoR experiment us¬ 
ing reasonable models of the field and the model beam of LOFAR 
produced by |Hamaker| ( |201 fj l usiim an electromagnetic simulation 
of the ASTRON Antenna Grounq and also test some leakage- 
correction strategies. The paper is organized as follows. Section]^ 
revisits the mathematical formalism of a radio interferometer and 
describes the DI errors and the LOFAR beam-related DD errors 
within the context of this formalism. Formalisms used for calibra¬ 
tion, imaging, flux conversion, RM synthesis and power spectrum 
analysis are also described briefly. In section we describe the 
pipeline and setup of the simulations of extragalactic point sources 
and present three different results: effect of DI errors and the ac¬ 
curacy of self-calibration in this case, effect of DD-errors and a 
possible DDE correction strategy, and finally errors due to self¬ 
calibration with incomplete sky models. Pipeline, setup and results 
of the simulation of Galactic foreground are presented in section 
12 where we show the results of rotation measure synthesis and 
power spectrum analysis, compare the power spectra of the leakage 
and the expected EoR signal, and test a potential leakage removal 
method. In sectionj^we give a summary of the paper, discuss some 
of the assumptions and limitations briefly and, finally, list the major 
conclusions of this paper. 


2 FORMALISM 

2.1 Mathematical model of a radio interferometer 

Here, we give an outline of the mathematical model of a radio inter¬ 
ferometer and refer the readers to jHamaker et al.| ( |1996[ >; |Smirnov| 
( |2011a[ > for a detail description. 

Consider a quasi-monochromatic electromagnetic wave prop¬ 
agating through space from a single point source. Using the Carte¬ 
sian coordinate system xyz where the signal propagates along 2 
direction, the signal, at a specific point in time (t) and space, can 
be described by the complex vector £{x, y, t) and transformations 
(e.g. contaminations) of this signal along its path can be represented 
by 2 X 2 Jones matrices. Assuming all such transformations to be 
linear, a cumulative Jones matrix (J) can be constructed from the 
products of the matrices. The signal detected by our telescope will 


^ M. J. Arts; http://www.astron.nl 
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be the intrinsic signal multiplied by this cumulative matrix, mathe- 
maticall5[^f ' = Jf. 

The electric field represented by this vector hits an antenna 
of our interferometer that has two feeds, each one sensitive to a 
specific polarization state of the vector in case of a perpendicularly 
incident electric field. Let us assume that the p and q feeds are 
sensitive to the x and y polarization states of the signal respectively. 
The feeds convert the respective electric fields into voltages and this 
conversion can be expressed as yet another Jones matrix yielding 


V = J'f' ^ 




( 1 ) 


Let us denote this antenna as a and assume that there is another 
antenna in our interferometer denoted by b. Voltages from each an¬ 
tenna are fed to a correlator that cross-correlates them to create 4 
pairwise correlations that can be written as a 2 x 2 matrix, known 
as the visibility matrix. 


= (VaVf) 


f i'^apVhp) 


(Vpp 

Vpq\ 

xi'^aqVfjp) 

i^aqVlq)] 

KVqp 

Vqq) 


which is related to the electric field correlations according to Eq.[^ 
i.e. 


V ab — Ja 


f (e^e*) 
\fyel) 


{eye*y)J 


(3) 


Here * denotes a complex conjugate, H the conjugate transpose or 
Hermitian conjugate and {) the time averages. Polarized waves are 
best described by Stokes parameters and their relation with the cor¬ 
relations of the electric field components, for a linear experiment, 
can be written as ([Hamaker et al.|1996^ 


f{exe*) (e,re*)\ ^ / 7 + Q U + iV'\^ 

\{eye:} {eye;}) \U - iV I-Q)-^ 

where B is the brightness matrix. Therefore, Eq. [^becomes 

V.i, = J.Bjf 


(4) 

(5) 


which contains all effects along the signal path in the form of Jones 
matrices. The effect fundamental to all interferometers is the phase 
difference between the measured voltages Va and V;,. To account 
for the phase delays in Eq.[^ consider the interferometer to be sit¬ 
uated in a Cartesian coordinate system represented by u, v, w and 
the antenna a to be located at the coordinates Ua = {ua, "fa, Wa)- 
The phase delay between the baselines a and b then becomes 


where Llab = Ua — Ub', l,m are the cosines of the right as¬ 
cension and declination of the source respectively; and n = 
y/1 — P — rri^. If we take out the phase delay scalar matrices (77- 
Jones) from J for both antennae and express them as a single scalar 
associated with the baseline, then Eq. [^becomes 

V<,i, = 3a'BKab3? = (7) 

where = B77ai) is called the coherency matrix as it represents 
the spatial coherence function ( |Clark|199^ of the electric field for 
this particular baseline. 

If, instead of a single source, we have a continuum of sources, 
the visibility matrix has to be written as an integration over all di¬ 
rections within the field of view and the cumulative Jones matrix 


® In this paper vectors are represented by calligraphy, matrices by bold and 
scalars by normal typefaces. 


has to be separated into two different matrices, one representing 
the direction independent effects (DIE, G-Jones) and another the 
direction dependent effects (DDE, 7?-Jones), 


^ab = Ga 


_Z,m 


Gf. 


( 8 ) 


This is the standard equation to describe the mathematical model 
of a radio interferometer that, from now on, we will refer to as the 
measurement equation. 


2.1.1 Mueller formalism 

For understanding the effects of systematic errors on the images 
produced from the visibilities, it helps to write this equation in 
terms of baseline-based Mueller matrices (M) instead of antenna- 
based Jones matrices (J) remembering the relation between the two 
( [Hamaker et al.|199^ , 

= jf)S (9) 


where the coordinate transformation matrix. 


S 


/I 1 0 0 \ 

1 0 0 1 i 

2 0 0 1 -i 

\1 -1 0 0 / 


( 10 ) 


and (g) denotes the Kronecker product. To do so, instead of tak¬ 
ing the matrix product of Va and v;, like in Eq. we have to 
take their Kronecker product to get the voltage correlation vector 
Vab ~ (Vpp Vpq Vqp Vqq)'^ where T represents transpose. Then 
Eq. [^becomes 

Vab = EabSlKab^^ dD 

l,m 


where Gat = Ga (g) G^, Eat = Ea (g) E^ and brightness vector 


2.7.2 Stokes visibilities 


In order to describe the relation between Stokes parameters and 
voltage correlations in Fourier space, let us define 




(at) _ 


= 3aZKab3b 


( 12 ) 


where = Vi,Vq,Vu,Vv is a Stokes visibility and Z = 

I, Q,U,V is a Stokes parameter. Comparing equations j 1 2| [7] and 
and remembering the definition of the coherency and brightness 
matrices, we can establish the relation between Stokes visibilities 
and the voltage correlations as (jSault et al.|1996[[Bunn|2007y 


dr = liVpp + Vqq) 

Vq = liVpp - Vqq) 

VU = liVpq + Vqp) 
VV = ^.{Vpq - Vqp) 


(13a) 

(13b) 

(13c) 

(13d) 
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Figure 1. Schematic diagram of a 24-tile LOFAR HBA station. A tile is 
made of 16 dual polarization dipoles. Dipoles see almost the whole sky 
(FWHM~ 90°), while the FWHM of a tile beam is ~ 20° and that of a 
station beam is only ~ 4°. There is a 15 cm gap between the tiles which is 
not shown here. 


2.2 Systematic effects 

In this section, we will discuss the effects of the systematic errors 
(G and E Jones) on the Stokes visibilities and the Stokes parame¬ 
ters for the case of LOFAR, although the aforementioned formal¬ 
ism is universal. LOFAR is a phased array covering the frequency 
range from 10-240 MHz. LOFAR stations consist of two types of 
antennae— LBA (low band antenna; 10-90 MHz) and HBA (high 
band antenna; 110-240 MHz). We use the HBA stations in our sim¬ 
ulations and a schematic diagram of a typical 24-tile LOFAR HBA 
core (situated within the central 3.5 km) station is shown in Fig.[^ 
In this case, 16 dipoles are combined to create a tile and 24 tiles 
are combined to create a station (for details see | van Haarlem et al.| 
|20T3l >. 


2.2.1 Direction independent effects 


To simplify calculations, while discussing DIEs, we will ignore the 
DDEs by assuming the i?-Jones terms of Eq.[^to be identity ma¬ 
trices. Consequently, the Mueller-matrix form of the measurement 
equation (Eq. [3 becomes, 


Vaf, = G, 


STKab^^=Gab 


SVz 


dldm 


(14) 


where Vz = 2LKab represents the Stokes visibilities without any 
systematic errors. The DIEs, denoted here by Gab, are caused by 
errors in the electronic gains of the antennae (gain errors) and non- 
orthogonal and/or rotated feeds (feed errors). Gain and feed errors, 
for antenna a, can be modelled by the Jones matrices, 

- ft' ft) ■“> = (-ft ‘7) <‘=> 


where gap is the gain error of the feed p of the antenna a and 
Cap is the spurious sensitivity of the p feed to the y polarization. 
The Jones matrix for all DIEs, i.e. G-Jones of Eq.[^ then becomes 
Ga = GaGft Gain and feed errors affect different Stokes visi¬ 
bilities (Eq. |13^ in different ways which can be illustrated by tak¬ 
ing into consideration how the Stokes visibilities observed by an 
instrument with DIEs differ from that of an error-free ideal instru¬ 
ment. Let’s assume that both G® and G^ of the ideal instrument 
are identity matrices and for a realistic instrument gains and feeds 
are in error by, 

= (ft" aL) ““ - (-ft ft) ■ 


Then, seven error parameters (hereafter Dl-error parameters) can 
be defined following|Sault et^jl996| equations 36-42) as. 


= (^pap ~\~ ^Qaq) ~\~ (^p6p ^Qbq) 

(17a) 

~ ~ ^9aq) + {^9bp ~ ^9bq) 

(17b) 

^U,V ~ i^9ap ^9aq) i^9bp ^9bq) 

(17c) 

^Q,U = (Cap + taq) + (Cfop + t^q) 

(17d) 

Sl,U = (Cap — ^aq) + {^bp ~ ^bq) 

(IVe) 

5l,V = (tap + taq) — (Cfop + Cfog) 

(17f) 

^Q,V ~ (tap tag-) {^bp ^bq) 

(IVg) 


where the subscript 7, Q stands for mixing between Stokes I and 
Q. Now, if the difference between the ideal Stokes visibilities and 
the Stokes visibilities affected by these errors is AV = — 

Vafj, then by assuming errors to be very small it can be shown that 


(see|Sault et al.|19961 appendix 

B), 




/ 5a 

Si,Q 

Si,u 

—i5i,v\ 


AV = -1 

Si.Q 

Si,u 

Ss 

—Sq,u 

^Q,u 

5s 

—iSq.v 

i5u,v 

> jideal 
^ab 


\—iSpv 

5q,v 

—i5u,v 

Ss y 



(18) 

Here, the 4x4 matrix is the instrumental Mueller matrix for the 
DIEs (hereafter DI-Mueller) and it determines the full Stokes re¬ 
sponse of an instrument without any direction dependent errors. 
It can be seen from the equation that a completely unpolarized 
source (Q, U,V = 0) will appear to have non-zero Stokes Q, U 
and V in an interferometric observation because of the DIEs Si,q, 
Sj^u and Sj^v respectively and these same errors will cause leak¬ 
age into Stokes I from Stokes Q, U and V respectively. The DIE- 
parameters can be used to determine calibration errors if, instead of 
comparing the ideal and the actual gains, we compare the input and 
the solved gains (|Sault et al.|1996|). 


2.2.2 Direction dependent effects 

Direction dependent errors in a radio interferometer are caused 
mainly by the Earth’s ionosphere and the primary beams—i.e. the 
radiation patterns—of the antennae. Here, we restrict ourself only 
to the LOFAR beam errors. The beam we use for the bowtie dipoles 
has been modelled by an analytic expression whose coefficients are 
determined by fitting to a numerically simulated beam raster gen¬ 
erated by the ASTRON Antenna Group ( |Hamaker|201 1[ hereafter 
HI 1). Here, we will give a brief overview of this model; for further 
details we refer the readers to HI 1. 

From basic symmetry considerations a generic expression for 
a dual dipole antenna i?-Jones matrix has been derived by Hll 
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Figure 2. (a) Direction dependent Mueller matrix representing the polarization response of the baseline 0-1 (127 m) of LOFAR at 150 MHz over the 3C196 
field (20° X 20°) at the time when the centre of the field culminates, (b) Spatio-temporal profiles as a percentage of total intensity—i.e. first row, first column 
(JV7]^]^) of the matrix representing Stokes I —for leakages from (1) I to linear polarization (7^), i.e. C^) linear to 7, i.e. + (3) / 

to circular, M 14 and (4) circular to I, M 41 . Here, A6 represents distance from the phase centre. See section[2.2.2|for details. 


which, for azimuth 0 and zenith angle 6 = (7r/2— elevation) can 
be written as, 

N 

Ee(e,<^) = y] R(fc',0)Pfc(0) (19) 

k'=0 


where the azimuth dependent rotation matrix 


R(fc',0) 


(cos[(—1)'“ {2k' + 1)<^] — sin[(—1)*’ {2k' + 

(sin[(-l)'''(2A:' + l)(^] cos[(-l)'''(2fc'+ 1)(^] j 

( 20 ) 


and the zenith angle and frequency (u) dependent projection matrix 
that contains the detailed geometry of the dipoles and the ground 
plane is 


'Py{e,v) 


(PB,k'{d,v) 

\ 0 


0 ) 

-P<t,,k'{0,v)) ’ 


( 21 ) 


and k' — 0 gives the ‘ideal’ beam, whereas the higher order 
terms represent the differences between the ideal and the more re¬ 
alistic beams. Each element of the projection matrix p{9, v), for 
each harmonic k', is calculated as 9[C\i> where 0 is a row vector 
(6»° ... 9^0), iz is a column vector {v° ... [C] is a 

2D matrix of dimensions {Ne -f 1) x {Nv + 1) that contains the 
complex coefficients determined by fitting to an electromagnetic 
simulation, and Ng = Ni, = 4. 

In Eq.[T^ Ee has been expressed in a topocentric (azimuth- 
zenith angle) coordinate, but in reality the source is carried around 


through the beam by the apparent rotation of the sky during 
an observation. To account for this effect, the position of the 
source is transformed from equatorial celestial coordinate system 
to the topocentric system. Eor polarized sources, there is an addi¬ 
tional factor— the relative rotation between the equatorial and the 
topocentric grids at the position of the source that causes the beam 
to rotate with the parallactic angle, known as the parallactic rota¬ 
tion which has been incorporated in the dipole beam model as a 
separate Jones matrix. Hereafter, by Ee we will refer to an element 
beam where all these effects have been taken into account. 

In an element beam Jones matrix the diagonal terms deter¬ 
mine the primary beam of the element and the off-diagonal terms 
the level of cross-polarization. Errors related to antenna pointing, 
beamwidth and beam ellipticity are all included in the diagonal 
terms. For a dipole of size D ~ 1.25 m the FWHM at 150 MHz be¬ 
comes A/D ~ 90° and the shape of the diagonal terms of the ma¬ 
trix is similar to an Airy pattern. The polarization response of a LO¬ 
FAR station is completely determined by Ee. Therefore, it would 
be interesting to analyse the beam Mueller matrix corresponding to 
an interferometer constructed by two such elements before entering 
into the discussion of the tile and the station beams. 

In a two-element interferometer, the component at the first row 
and first column of the Mueller matrix (hereafter Mn) represents 
the Stokes / response of the interferometer to a completely unpo¬ 
larized point source of unity flux and M 12 gives the corresponding 
Stokes Q response. Examples of Stokes I and Q responses of a 
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Figure 3. (a) Gaussian fit to the azimuthally averaged Stokes I response of the 0-1 baseline of LOFAR at 150 MHz over the 3C196 field when the field 
culminates (Mu component of Fig. |^). (b) FWHM of the Stokes I beam at different frequencies (solid); the a\/D curve (dashed) is overplotted, (c) A 
single line through the centre of the Mueller term responsible for linear polarization leakage (see caption of Fig.|^ at different frequencies. The leakage is 
shown as a percentage of Stokes I flux density. 


LOFAR LBA dipole can be seen in Fig. 3.8 and 3.9 of |Bregman| 
( |2012[ hereafter B12) respectively. From the figures we see that 
Stokes I response is almost circular with amplitudes decreasing 
from the centre toward the edges until the first null. Stokes Q re¬ 
sponse, on the other hand, has a cloverleaf pattern with 2-fold sym¬ 
metry corresponding to the physical structure of the dual dipole. 
The cross-polarization over a heam is conventionally measured hy 
the ratios Q{0)/I{0), U{6)/I{6) and V{9)/I{0). Comparing Fig. 
3.9 and 3.8, B12 finds that Q// is lowest at the centre and increases 
quadratically with 9 and reaches a value of 0.5 at the FWHM. It 
implies that an unpolarized source situated at FWHM of a dipole 
heam will become 50% polarized in the observed data due to in¬ 
strumental polarization. 

The beams of the 16 dipoles (Ee) in a tile are combined in an 
analogue way to form the tile beam which is narrower (~ 20°, Fig. 
[T]( and the beams of all the tiles in a station are digitally combined 
to form the station beam which has the smallest width (~ 4°). As¬ 
suming the tile beams (Et) have been created by phasing the con¬ 
stituting dipole beams, the beam of the station a can be written as 
( |Yatawatta|2009^ 

E4e,0) = w"v(k)©Et(6l,0) (22) 


where © denotes the Hadamard product, k is the wave vector, v(k) 
is the steering vector, i.e. the delay an incoming wavefront expe¬ 
riences depending on the position (rt) of the observing tile in a 
station that can be expressed as 


r(k) = 




o-fk.rjv-1 


(23) 


for N number of tiles and w is the weight vector that contains the 
complex weights associated with each tile. Station beams cut only 
a small portion of the element beam and get a polarization response 
depending on which part of the element beam it is tracing. The side- 
lobes of the station beam cut yet another part of the element beam 
and accordingly acquire a different polarization response. Station 
beams that are formed to track a source in the sky follow a trace in 
azimuth and elevation over the polarized element beam. Hereafter, 
by beam we will refer to the beam of a single station, Ea. 


We could, in principle, derive a direction dependent equiva¬ 
lent of Eq. [fusing Ea as the only systematic error and ignoring 
the DIEs, but it will be much more complicated in this case. So, 
instead, we numerically calculate the baseline-dependent Mueller 
matrices (e.g. Eab) from the constituent station beams (Ea and Eb) 
following the formalism of section [?.l.l[ Such a Mueller matrix for 
baseline 0-1 (a 127 m baseline formed by the two sub-stations of 
the central core stations, CSOOIHBAO and CSOOIHBAI) at 150 
MHz, at the time when the centre of the target field (20° x 20°) 
culminates has been shown in Eig.|^. The components of the ma¬ 
trix have been normalized with respect to the Mueller matrix at the 
phase centre resulting in a differential Mueller matrix; hereafter, 
by differential beam or nominal beam we will refer to this form 
of the Mueller matrix. Let’s denote this matrix by where the 
superscript represents the station numbers. 

can be thought of as a direction dependent equivalent 
of the DI-Mueller (Eq. [T8l >, hence we can call it the DD-Mueller. 
By comparing these two matrices, we can see that M 21 component 
of the DD-Mueller will cause Stokes I to leak into Stokes Q. The 
off-diagonal terms of show the spatial variation of the instru¬ 
mental polarization— it is lowest at the phase centre and increases 
toward the edges until the first null and then, after a gap, we get 
further polarization at the location of the first sidelobe. In addition 
to the spatial variation, all components of the instrumental Mueller 
matrix also vary with zenith angle, or equivalently with hour angle, 
of the source during an observation. To show the dependence on 
the directions and sidereal time simultaneously, i.e. spatio-temporal 
dependence, we calculated for all hour angles. In Eig.l^, we 
show spatio-temporal profiles of various leakages as a percentage 
of total intensity. Leakage from linear polarization to total inten¬ 
sity, i.e. \/ATj 2 + X 100, at different distances from the 

phase centre (x axis) and at different hour angles (y axis) during an 
eight-hour observation is shown in the top panel. The second panel 
shows fractional leakage from Stokes 1 to linear polarization and 
the third and fourth panels show fractional I ^ V and H —> / 
leakages respectively. These figures show the variation of the leak¬ 
ages along a single line through the centre of the field at every hour 
angle during a night-long observation. 

From the spatio-temporal profiles, we see that leakage in¬ 
creases with both distance from the phase centre and zenith angle. 
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During the beginning and the end of the observation zenith angle 
is very high and the beam is extremely attenuated which results in 
a very high percentage of leakage. Leakages vary across the FoV 
mainly due to polarization aberrations caused by geometric projec¬ 
tion of the antenna on the plane perpendicular to the line of sight 
(see section 5.3 and Fig. 2 of |Carozzi & Woan|2009^ . The projec¬ 
tion changes as a function of direction and zenith angle because of 
both the coordinate rotation and parallactic rotation that were intro¬ 
duced in the beam model (as discussed before). We see that at high 
zenith angle the leakages change more rapidly, but these effects can 
be considered constant within ten minutes (B12) which is an useful 
assumption for primary beam correction. 

Besides direction and elevation, the width and shape of the 
beam also vary with frequency. Fig.[^ shows a Gaussian fit to the 
azimuthally averaged station beam (Mu) that gives us an FWFIM 
of 3.8° at 150 MHz. Fig. shows the beamwidths obtained by 
Gaussian fitting as a function of frequency and we can see that the 
curve closely follows the a\/D relation where A and D denote 
wavelength and station size respectively (for an analogous fitting, 
see Fig. 21 of |van Haarlem et al.|2013^ . Leakages also vary with fre¬ 
quency, albeit not in a very prominent way; as evident from Fig.[^, 
within approximately ten degrees leakage changes very slowly with 
frequency. Therefore, if we have multi-frequency data, the leakages 
can be removed by utilizing their spectral smoothness. 

Ideally, the beam should be exactly same for all elements 
and, consequently, for all baselines, for traditional calibration to 
work efficiently, but making them slightly different in configuration 
could be advantageous in another way. In case of LOFAR, although 
all dipoles are rotated into the same position, station configurations 
are rotated with respect to one another to minimize blind angle ef¬ 
fects and to average out the effect of grating lobes (B12). 

2.3 Calibration and imaging 

In Dl-calibration, it is assumed that all baselines of an array observe 
the Fourier transform of a common sky which is only true if DDEs 
are taken to be identical across all antennae. Consequently, of 
Eq. [^becomes a function of just I, m and the common sky observed 
by all baselines becomes = EBE^, i.e. the true sky attenuated 
by the beam. Then, Eq.j^can be written as 

= G.X°i,Gf (24) 

where X°i, is the element by element 2D Fourier transform of Be. 
The most widely used Dl-calibration method, self-calibration or 
selfcal works with this form of the measurement equation. The 
first step of selfcal is to create a model of the observed sky and 
to ‘predict’ the corresponding visibilities, that an interfer¬ 

ometer would produce. Then, the values of G terms that minimize 
— Vai, are determined. G terms can be calculated to a very 
high accuracy, because an array provides over-determined informa¬ 
tion as N{N — 1) complex visibilities are available for computing 
only 2N — 2 error parameters, N being the number of antennae. 

The inferred values (G) are applied to the observed visibilities 
to yield the corrected visibilities as 

V:r = G-^V^tG,-". (25) 

Inverse Fourier transform of the weighted and gridded visibilities 
produce a ‘dirty’ image, which is the true sky convolved with the 
PSF. To recover the true sky as sampled by the visibilities as closely 
as possible, the PSF is deconvolved from the dirty image iteratively 
producing a ‘clean’ image. As the primary beam has not been cor¬ 
rected for, this clean image is actually the true sky attenuated by 


the primary beam (Be). If the primary beam is assumed to be same 
for all antennae and at all times, the true brightness distribution B 
can be extracted from Be by just multiplying it with the inverse of 
E. Traditionally, this is what has been done for dish instruments 
with small FoV. But in case of wide FoV instruments, e.g. LOFAR, 
time-frequency-baseline variations of the instrumental Mueller ma¬ 
trices (M, Fig. cannot be ignored and one way of dealing with 
this is AW-projection ([Tasse et al.|2013|>. 


2.3.1 AW-projection 


The problem of imaging can be expressed in Mueller formalism as 
V = Al + e where V is the total set of visibilities, I is the set of 
Stokes images to be estimated, e is the noise, A — WSJ-A4 ig¬ 
noring the ionospheric effects, W is the set of visibility weights, S 
is the sampling function, J^ is the Fourier transform kernel, and A4 
is the Mueller matrix corresponding to the primary beam. Each of 
these parameters is a multi-dimensional matrix (for explanation see 
|Tasse et al.|2013^ . AW-projection, as implemented in AWImager, 
calculates I, an estimate of I, iteratively as, 

i"+i =I"-f ( 26 ) 


where $ is a non-linear operator that estimates the deconvolved sky 
from the residual dirty image (V — ^I"). Here the construction 
of the residual dirty image constitutes the major cycle and the de- 
convolution the minor cycle. Note that .41" is the forward Fourier 
transform taking into account all instrumental effects and this has 
to be done accurately for the solutions to converge; during predic¬ 
tion of visibilities using AWImager, only this step is performed. 
On the other hand, during minor cycle only an approximation of 
{A^A)~^ is calculated and applied on the residual. A-projection, 
as described in [Bhatnagar et al.|j2008| l, is a fast way for applying A 
or A^. In AWImager, the element beam (Ee) and the array fac¬ 
tor (w^v(fc) of Eq. 22 1 of LOFAR have been taken out of the A4 


matrix of the A-term and they are applied separately. 


2.4 Flux conversion 


For easier comparison with the predicted level of the EoR signal 
we convert fluxes to intensities and express them as temperature. 
If Fjy is the flux of a radio source in Jy, then the corresponding 
intensity in K units can be written as. 


Tk 


^^Pjy ]^q-26 
2hB^E 


(27) 


where ks is the Boltzmann constant and Q,e = / (4 In 2) is the 

beam solid angle, 9 being the FWHM of the Gaussian restored PSF 
calculated during imaging. 


2.5 Rotation measure synthesis 

The rotation of the plane of polarization (x) of a linearly polarized 
signal while propagating through a magnetized plasma is called 
Faraday rotation which, for a single Faraday screen along the LOS, 
can be written mathematically as x = Xo + where xo is the 
intrinsic polarization angle and Faraday depth, 

^observer 

<E> = 0.8l/ neB||d( (28) 

J source 

where ris is the density of electrons and By is the magnetic field 
component along the LOS. Note that rotation measure (RM) is de¬ 
fined as dx/dJ? and hence for a single phase screen along the LOS 
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it is equivalent to Faraday depth. Polarized surface brightness per 
unit Faraday depth, F{^) can be obtained from the polarized sur¬ 
face brightness per unit squared-wavelength, f’(A^) using the tech¬ 
nique of RM-synthesis jBrentjens & de Bruyn|2003| l; mathemati¬ 
cally, 

/ °° 2 

d\^ (29) 

- OO 

where is the Fourier transform of the wavelength sampling 

function, known as ‘rotation measure spread function’ (RMSF) and 
* denotes convolution. 

The polarized brightness, V = Q + iCij^is a complex valued 
function and, hence, F{^) is also complex. However, a Faraday 
dispersion function for real valued Stokes I, F/(‘f>) can also be 
calculated assuming its imaginary parts to be zero in all spectral 
bands (e.g. |Geil et J^|2011) . As the Fourier transform of a real 
function is always Hermitian, Fj ($) = Tj(—<&). The same can 
be done for Stokes V. In section]^ we will present some of our 
results in terms of T’($), TV(3>) and TV('I?). 


2.6 Power spectrum analysis 

The power spectrum (hereafter PS) of an image is the measure of 
the variance per unit angular wavenumber (k — 2 -k/6). As the first 
detections of the EoR signal will be statistical, and its PS is the most 
widely used statistic (e.g. JBowma^t alj2006 yHarker^t^.|20 1 0[ 

[Moore et al.|2013] |Patil et al.||2014[ |Chapman et al.|2014^ , most 

of our analysis will be done through PS. We present three types 
of PS: 2D, 3D cylindrical and 3D spherical, and in all of them the 
wavenumbers are converted to the unit of comoving Mpc”'^ at the 
redshift corresponding to the observing frequency. PS can be calcu¬ 
lated from the weighted visibilities directly. As the imaging process 
puts weights on the visibilities and calculates the resulting PSF, we 
have measured the PS from the Fourier transform (hereafter FT) of 
the images remembering that the squared complex modulus of a FT 
yields the PS of a signal. 


2.6.1 2D power spectrum 


Assume that is the 2D FT of the image Iim where u, v rep¬ 
resent the spatial frequencies corresponding to the angular scales 
I, m. The minimum and maximum spatial frequencies of luv are 
determined by 1 /{N^OpiA and l/(20pix) respectively where Spix 
is the angular size of the pixels in Iim and 
where Ni and Nm are the total number of pixels in I and m di¬ 
rections respectively. We cut the portion of delimited by the 
minimum and maximum physical baselines and calculate the 2D 
PS asP 2 D(M,u) = 

To produce ID angular PS, we divide P 2 D in several concen¬ 
tric circular bins and calculate the average power at every bin. Fi¬ 
nally, we plot the average power in the bins as a function of co¬ 
moving transverse wavenumbers corresponding to the bins defined 
as ([Morales & Hewitt|2004| equations 2-3), 


fc_L 


2ttUx 

DAz) 


(30) 


where Ux = in units of wavelengths, transverse comov¬ 

ing distance at redshift «, Dc{z) = dz'/E{z'), and dimension¬ 
less Hubble parameter, E{z) = [Hm(l + 2 )^ + Hfe(l + z)^ + 
Dm, Hfe and Ha being the matter density, curvature and 
cosmological constant parameters respectively. Thus, we obtain k± 
in units of Mpc“^ and P 2 D{k±) in units of Mpc^. Note that the 
minimum and maximum values of k± are determined by and 
respectively, as shown in[Vedantham et al.[([2012[ equations 
13-14). 


2.6.2 3D power spectrum 

Assume that luvr/ is the 3D FT of the image Iimv where rj rep¬ 
resents the LOS spatial frequency corresponding to the LOS dis¬ 
tance signified by the frequency 12 (see [Morales & Hewitt[|2004| 
Fig. 2). After taking only the portion of the cube that represents 
real baseline distribution as before, the 3D PS can be calculated 
as P3d{u, V, rj) = Two types of binned PS can be calcu¬ 

lated from this PS-cube; cylindrical, P 3 D{k±,k\\)^ spherical, 

P3D{k). 

In the cylindrical case, averaging is done in concentric cylin¬ 
drical bins centred on the centre of the cube. Hence, P 3 D{k±, A:||) 
is the average power of all uv cells within a logarithmic cylindrical 
bin around k±, fcy where the comoving LOS wavenumber. 


fell =77 


2-kHoE{z)u2i 
c(l -I- 2)2 


(31) 


1221 being the rest frequency of 21-cm radiation emitted by HI, and 
k± is the same as defined by Eq.|^ The minimum and maximum 
values of fey are given by tymin = 1/P and 77max = 1 /Ai 2 respec¬ 
tively where B is the bandwidth and Au is the frequency resolution 
provided by the instrument. From the minimum and maximum val¬ 
ues of k± and fey, it is evident that the boundaries of the fc-space 
are defined by the instrumental parameters (see e.g. |Vedantham etj 
[ar][2012| Fig. 4). Instead of showing the raw power we plot the 
quantity A‘^{k±,kA = k\k^^P3D{kj_, k\\) /in our 2D fig¬ 
ures which has the dimensions of temperature squared. 

For constructing the spherical 3D PS, we divide the PS-cube 
in concentric spherical annuli around the centre of the cube and 
average the power in every annulus. Consequently, we get a ID 
PS as a function of fc = Pi°i quan¬ 

tity Ad{k) = k^P 3 o{k) / {2'k‘^) that has the same dimensions as 
A"(fcA,fcy). 


3 SIMULATIONS OF EXTRAGALACTIC FOREGROUND 

To show the effects of direction independent errors on calibration, 
we simulate the observations of a mock sky with point sources. In 
case of the direction dependent errors, we first simulate a mock sky 
to show the trend of the effects, and then proceed to simulate the 
realistic sky to quantify the effects expected in the LOFAR-EoR ob¬ 
servations. We did not include any additive noise in the simulations 
described in this section. Below we describe the general pipeline of 
the simulations followed by the set-ups and results of the specific 
simulations. 


3.1 Pipeline 

® In this paper P always refers to Q + iU, while P is always \Q + iU\, 

and note that the 2D and 3D power spectra, denoted by P 2 D and P 30 A block diagram of the pipeline for simulating extragalactic point 

respectively, are not related to P or P. sources is shown in Fig. We start from a given model of the 
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Figure 4. Block diagram of the pipeline of the simulations of extragalactic 
foreground. Blocks with solid and dotted borders represent simulations with 
DD and DI errors respectively; blocks with dashed borders represent steps 
performed for both simulations, but separately. Airows with dashed line- 
styles have been used to avoid intersection between arrows. FT, IFT and SC 
stand for Fourier transform, inverse Fourier transform and self-calibration 
respectively. 

sky (described in the specific sections) and predict the visibilities 
that LOFAR would produce in the presence of certain DI and DD 
(beam) errors. Simulations with the two systematic errors are done 
separately, although some steps are common to both of them. 

DI errors are introduced in accordance with the formulation 
described in section |2.2.1| After prediction, the visibilities cor¬ 
rupted by the DIEs are self-calibrated using the same sky model 
that was used to predict. Then, the gains determined by selfcal 
are compared with the input gains to calculate the error parame¬ 
ters defined by equations |17^-g. Additionally, the solved gains are 
applied to the model visibilities to produce corrected visibilities. 
All processes up to this point are performed using the standard 
LOFAR calibration and simulation software. Black Board Selfcal 
(BBS; [Pandey et al.|2009[ l. We image both the corrupted and the 
corrected visibilities using CASA and produce 2D PS from the im¬ 
ages through the procedure described in section [2.6.1| 

DD errors are introduced by multiplying every point source 
in the model with the relevant station beam at the position of the 
source at every timeslot. Fourier transform of the beam attenuated 
sky yields the visibilities corrupted by DDEs. We carry out two 
different simulations with these dataset: one to measure effects of 
DD errors, and another to quantify the errors in calibration due to 
incomplete calibration sky model. The latter could be done mean¬ 
ingfully without introducing systematic errors at all, but we did it 
this way to make it more realistic. 

To quantify the effects of DD errors, first, we correct the cor¬ 
rupted visibilities for the beam at the phase centre which, in real¬ 
ity, normalizes the DDEs with respect to the phase centre so that 
only the differential nominal beam effects remain (this step is not 


Table 1. Observational setup for simulations of extragalactic sources: 


Number of LOFAR HBA stations used, N 

59 

Number of baselines, N{N — l)/2 

1711 

Number of spectral subbands 

1 

Number of channels in the subband 

1 

Central frequency of the channel 

150 MHz 

Width of the channel, i.e. frequency resolution 

0.19 MHz 

Total observation time 

8 h 

Integration time, i.e. time resolution 

10 s 

Number of timeslots 

2874 

Number of visibilities 

5090520 

Baseline cut (Wmin ~ ■^^max) for imaging 

0.06 - 20 km 

Baseline cut for PS estimation 

0.06 - I km 

Angular resolution (PSF) of the images, aA/iimax 

~ 0.34 arcmin 

Physical width of the HBA stations, D 

30 m 

FWHM of station primary beams, aXfD 

~ 3.78 deg 

Field of view, 7r(FWHM/2)^ 

11 .2deg2 


shown in Fig. [^. Then, we image both the corrupted and uncor¬ 
rupted (ideal) visibilities and produce 2D PS from the images. Fur¬ 
thermore, we extract the fluxes and positions of the brightest point 
sources in the corrupted and uncorrupted images using PyBDS]v(3 
and compare them. Finally, we correct the visibilities for the dif¬ 
ferential beam and produce images from them using AWImager. 
Fluxes of the beam-corrected images are compared with the uncor¬ 
rected fluxes to quantify the quality of the correction. 

To determine calibration errors due to an incomplete sky 
model, we calibrate the corrupted visibilities using different incom¬ 
plete sky models. As the same DDEs are included during both pre¬ 
diction and calibration, the remaining errors will be only due to 
the incompleteness of the models. The deviation of the different 
corrected visibilities from the corrupted visibilities is demonstrated 
through PS. 

3.2 Direction independent errors 

To show the effects of DI errors and test their correction strat¬ 
egy, we ignore the DDEs and introduce DIEs for every station and 
timeslot as G-Jones matrices. Both gain (g) and feed (e) error terms 
of G are modelled as complex numbers that are random at every 
time-step drawn from a Gaussian distribution with zero mean and a 
certain standard deviation (rms). Then, we create a sky model con¬ 
taining 25 sources of 5 Jy Stokes / flux (Q, U,V = 0) in a 5 x 5 
uniform grid of 1° separation, predict the DIE-corrupted visibili¬ 
ties for all baselines of LOFAR and perform all the other steps de¬ 
scribed in the previous section and shown in Fig.|^(see the blocks 
with dotted and dashed borders). The rms of the introduced errors is 
the same for every term of the G-Jones of every station and we re¬ 
peat this experiment thrice for three different rms Dl-errors; 10“^, 
0.01 and 0.1. Note that, as the calibration was done with a perfect 
sky model, the errors will be due only to the calibration process 
itself. 

We analyse the results using two parameters: fractional rms 
selfcal error (5f ) and square-root of the residual power spectrum 
{\/P{k)) which, in effect, gives the rms of the images at different 
spatial frequencies. To determine 5f, we calculate AG® and AG-^ 
(see Eq. [T^ by differencing the model gains and the solved gains 
for two stations, and then, calculate the DIE-parameters (5) for the 

http ; //t inyurl. com/PyBDSM-doc 
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Figure 5. Left: Fractional error on the 7 DIE parameters defined in Eq. |17| for a single baseline as a percentage of the input rms Dl-en'ors. The Ag and Ae 
used to calculate these parameters are the differences between the components of the input Jones matrix and the Jones matrix calculated by self-calibration. 
Right: Square-root of the fractional residual power spectra, which is equivalent to the rms of the image, for different rms Dl-errors and Stokes parameters. See 
section lr^ for details. 


baseline created by those stations (equations |17t. We did not plug in 
the values of AG® and AG^ directly in Eq.|l6|to calculate 5, but 
created an error DI-Mueller matrix from the AG matrices of two 
stations following Eq. and extracted the 7 relevant parameters 
from it. 5f for a given 5 is the rms of the <5 as a percentage of 
the input rms Dl-error. The seven 5f are plotted as a function of 
the input rms Dl-errors on the left panel of Fig. We see that 
fractional selfcal errors increase linearly with rms Dl-errors, and 
for an rms Dl-error of 10“^, which is not unrealistic, the error on 
these parameters is less than 0.002%. 

For calculating residual P{k), we subtract the corrected 
Stokes images from the corrupted ones and measure the PS of the 
residuals. As we did not subtract any source from the corrected vis¬ 
ibilities, if the calibration error is low the difference between the 
corrected and the corrupted visibilities should also be low. Pr{k) is 
the PS of a residual image as a percentage of the PS of a corrupted 
image, of the different Stokes images for the three simu¬ 

lations are plotted on the right panel of Fig. which clearly shows 
that the calibration errors propagated to the PS are negligible as 
expected in the absence of additive noise. For an rms Dl-error of 
10“^, errors on y/F^{k) or, equivalently, on the rms of the image is 
less than 0.005%. Furthermore, by comparing the Stokes I,Q + iU 
and V power spectra for an rms Dl-error of 0.01, we see that the er¬ 
rors on different Stokes parameters are the same, as expected. This 
simulation shows that self-calibration can correct for the Dl-errors 
to a very high accuracy if we have a sufficiently accurate model of 
the sky. 


3.3 Direction dependent errors 

To show the effects of DD errors on point sources and to test one 
of their correction strategies, we ignore the DIFs, introduce DDEs 
as station beams and carry out the steps outlined in Fig. (see 
the blocks with solid and dashed borders). As mentioned before, 
we implemented two different simulations with the DDE-corrupted 
dataset; the purpose of the first one is to show the effects of DD- 
errors on the Stokes parameters and this has been done for two 
different sky models, some information about which are listed in 
table 12 


3.3.1 Test with a mock sky 

To show the general trend of the effects of DD-errors, we make 
a mock sky model comprising 225 unpolarized point sources ar¬ 
ranged in a 15 X 15 uniform grid of 0.66° separation centred on 
the position of 3C19tp]and simulate an 8-hour, 150-MHz observa¬ 
tion of LOFAR, taking into account the beams described in section 
|2.2.2[ The source at the centre of the grid is given a flux density of 
100 Jy, while each of the other sources have a flux density of 0.4 Jy. 
The central source has been made exceptionally bright (analogous 
to the 3C196 field) to be able to check the consequence of calibrat¬ 
ing an otherwise dim sky with a very bright point source which will 
be described in section [T4l 

As the sources were completely unpolarized, the Stokes 
Q,U,V images created from this dataset contain only the flux 
leaked from Stokes I, i.e. instrumentally polarized sources. These 
sources are shown on the middle and right panels of Fig. Each 
bubble in the plots represent an instrumentally polarized point 
source and the size and colour of the bubble represent the flux of 
the source as a percentage of its Stokes I flux. The figures show that 
leakages to both linear and circular polarizations increase as we go 
out from the centre of the field. As for the levels of leakage, within 
the central 4 degrees, i.e. within the first null of the primary beam 
at 150 MHz, linear polarization leakage (/ —>■ P) is around 0.5%, 
and circular polarization leakage (/ —>■ V) is less than 0.003%. In¬ 
strumental polarization of the central bright source (not shown in 
the figure) is very low, because before imaging the visibilities cor¬ 
rupted by the DDEs were corrected for the element beam (Ee of 
Eq.|19|l at the phase centre, thereby making the leakage terms very 
close to zero at that point. In physical terms this means that the pro¬ 
jection of the beams on the sky had been made perfectly orthogonal 
at the phase centre. What is left after this centre-correction is the ef¬ 
fect of the differential beam (e.g. Fig.|^). There is an anomaly in 
the south-east comer of the middle and the right panels of Fig. 
which can be attributed to the errors in extracting fluxes of very 
dim sources situated near the null of the primary beam. 

These results are consistent with the beam model described 


A quasar situated at z ~ 0.871 with a flux density of 74.3 Jy at 174 
MHz. 
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Position error [arcsec] /->Pieakage [%] /-^Kieakage [%] 
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Figure 6 . Left'. Distribution of 103 sources from the 225 sources arranged in a 15 X 15 uniform grid within a 10° field of view. Flux (bubble size) and position 
(colour) errors due to calibration with only the prominent central source are shown. Middle'. Same distribution with coiresponding fluxes leaked from Stokes I 
to linear polarization as a percentage of Stokes I flux (size and colour). Right: Same as the middle figure except that it is for the leakages to circular polarization 
which is much lower. 


Table 2. Sky models used for the different simulations of extragalactic foreground with DDEs. 


Field 

Phase centre 

Phase centre 

FoV 

Catalogue 

Number of 

Maximum 

Minimum 

Totapl 

Spectral index 


(Equatorial J2000) 

(Galactic) 

(deg) 


sources 

(Jy) 

(Jy) 

(Jy)^ 


Mock 

a ~ 8'*13'36", (5 ~ 48°13'0" 

1 ~ 171°, fc ~ 33° 

10 


225 

100 

0.4 

189.6 

-0.75 

3C196 

a ~ 8'*13'36", (5 ~ 48°13'0" 

1 ~ 171°, fc ~ 33° 

10 

nRSi]^ 

4567 

83 

0.027 

796.64 

-0.75 


“ All flux densities shown here are at 150 MHz. 

^ The Faint Images of the Radio Sky at Twenty-cm survey, produced by NRAO VLA at 1365 and 1435 MHz and 5" resolution; noise ~ 0.15 mJy. 


in section |2.2.2| For example, we can understand both the trend 
and the level of linear leakage seen in Fig. by comparing it to 
the M 21 and M 31 components of the instrumental Mueller matrix 
shown Fig. |^, or to the spatio-temporal profiles of the leakages 
shown in Fig.|^. We expect to see leakage at this level also in the 
realistic simulations and this expectation will be put to the test in 
the next section where we describe the simulation of one of the 
LOFAR-EoR target fields. 


3.3.2 3C196 field 

The 3C196 field (centred on the bright quasar, 3C196; [Bemardi et| 
|al.|2010^ is well-suited for EoR observations because the presence 
of a bright and almost unresolved source at its centre allows very 
accurate direction independent calibration, and it is situated in one 
of the colder regions of the Galactic halo. To make an unpolarized 
sky model for simulating this field, we extract Stokes I fluxes and 
positions of the sources brighter than 25 mJy within a radius of 5° 
around 3C196 from the FIRST survey catalogue (see table and 
extrapolate the fluxes to that of 150 MHz using a spectral index of - 
0.75 which is typical for the radio sources at these frequencies. The 
eponymous source, 3C196, has been taken out of this model, and 


a 4-component improved model of the source made from LOFAR 
data by V. N. Pandey has been inserted in its place. 

Linear leakage of the brightest 33 sources (Stokes I > 100 
mJy) is shown on the right panel of Fig. Both colour and size 
of the bubbles in the figure represent the percentage of leakage. 
Extraction of fluxes and positions of the sources in this case is not 
as precise as that of the gridded sky model as here sources are much 
more closely spaced; thus some errors in this scatter plot originate 
from the source extraction process. Nevertheless, the figure, as a 
whole, is quite informative; we see that linear leakage can be as 
high as 4%, but for most of the sources it is less than 2% and for 
the sources very close to the phase centre only less than a percent 
leak, as expected. The sources with the highest leakages (the three 
reddest bubbles) are very dim in Stokes I which can be seen by 
comparing these three bubbles with the corresponding bubbles on 
the left panel where the Stokes / fluxes are shown as colour of the 
bubbles. These leakages might not be real, but a consequence of 
errors in the source extraction process. The leakage from 3C196 
itself is very low and hence is not shown here. 

The overall level of the leakage can be better understood from 
the fractional (as a percentage of Stokes I PS) power spectra of 
Stokes Q, U, V shown in Pig.[^. The PS of Q/I and U/I tell us 
that the rms of the linear leakage is 0.05 ~ 0.06% of the rms of 
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Figure 7. Left. Distribution of the brightest 33 sources {I > 100 mJy) in the 3C196 field with their corresponding Stokes I fluxes (colour) and flux en'ors 
(bubble size) due to calibration with different number of sources (numbers in the legend) in the sky model. The percentages in the legend refer to the minimum 
and maximum flux errors. Note that after calibration with 1000 sources errors for most of the sources decrease. Right. Same distribution with the corresponding 
linear polarization leakages as a percentage of Stokes I flux. Both colour and size of the bubbles represent fractional leakage. 


the Stokes / image. On the other hand, rms of circular leakage is 
almost 4 orders of magnitude lower. Leakage from linear polariza¬ 
tion to Stokes I, which is relevant for the EoR experiments, will 
be similar to the I ^ P leakage shown in this simulation, as ev¬ 
ident from a comparison of the first and second panels from the 
top of Fig. w However, compact radio sources are usually un¬ 
polarized or very weakly polarized and hence the leakage from 
polarized point sources into Stokes I is very low and even that 
leakage can be removed by direction dependent calibration (e.g. 
SAGECal; |Kazemi et al.|2011[|Kazemi & Yatawatta|2013| > and/or 
AW-projection jTasse et al.|2013^ . 

3.3.3 Correcting polarization leakage of point sources 

There are several strategies for correcting beam-related DDEs 
which are classified broadly into two categories: image-plane and 
Fourier-plane corrections. Here, we test one of the Fourier-plane 
strategies called AW-projection (see section [2!3.1[ l, a particular ver¬ 
sion of which is implemented by AWImager for the LOFAR AW- 
terms. 

To do a simple test, we create a dataset from a sky model 
consisting of 36 unpolarized 10 Jy sources in a 6 x 6 uniform 
grid of 0.5 degree separation so that all sources are within the 
FWHM of the primary beam, and then try to correct the Stokes 
1 fluxes and remove the leakages using AWImager. As mentioned 
in section [2.3.1| LOFAR ^-terms are separated into two parts by 
AWImager: the slowly varying (in time) element beam (Ee), and 
the fast-varying array factor. We assume Eg to be constant within 
12 minutes, and the array factor to be constant within 5 minutes. 

The result is shown in Fig.j^ both color and size of the bubbles 


represent percentage of leakage removed by AWImager. It seems 
that up to 80% of the leakage can be removed. The performance 
appears to be worse near the centre of the field than further away 
which is counter-intuitive, but the leakage is already very low near 
the centre and the bad performance could be due to the inefficiency 
of both AWImager and the flux extraction software in dealing with 
faint sources. We should be careful to draw any final conclusions 
on the effectiveness of AWImager in removing leakages from our 
data as the software is still under construction and we are not aware 
of any test of leakage removal done on a realistic dataset. 


3.4 Selfcal errors due to incomplete sky model 

Incomplete sky models can lead to many problems in directionally 
independent self-calibrated data, among them generation of spuri¬ 
ous source components, removal of real source components and the 
generation of ghost sources, a spurious source whose flux is propor¬ 
tional to the flux of an unmodelled source ( [Grobler et al.|2014| . We 
try to quantify the calibration errors due to incomplete sky models 
with different numbers of sources in the models for the 3C196 field. 
Note that, as we included the beam during both prediction and cal¬ 
ibration, its effect was taken out and we were left with only selfcal 
errors. The calibration is performed using BBS which is based on 
the matrix formalism described in sectionj^ 

We make 5 different calibration sky models that contain 
roughly 10, 15, 30, 60 and 75% of the total flux of the field; the 
models have 1, 5, 50, 400 and 1000 sources respectively. After self¬ 
calibrating the field with each one of these models, we calculate the 
difference of fluxes of the sources between the calibrated and un¬ 
calibrated data, and also create corresponding residual PS. The left 
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Figure 8. a\ Square-root of the power spectra of the Q, U, V leakages as a percentage of the Stokes I power spectrum within the central 10 degrees of the 
3C196 field, b'. Residual (after subtracting calibrated data from the uncalibrated ones) power spectra of Stokes / as a percentage of the uncalibrated Stokes I 
PS of the same field. The different cases are for calibration with different number of sources in the sky model. These residuals correspond to calibration en'ors 
due to incomplete sky model. 


panels of Fig. show the flux errors on the sources that contribute 
to the largest errors in the field; the filled bubbles represent the 
errors after calibrating with only 10% of the total flux, while the 
unfilled bubbles with red borders are for the case when 75% flux 
is modelled. According to the figure, errors go down significantly 
after improving the sky models. 

In Fig.[^, we show the PS of the Stokes I residual after sub¬ 
tracting the calibrated images from the uncalibtrated ones as a frac¬ 
tion of uncalibrated Stokes / PS. As there is an exceptionally bright 
source (the second brightest source is only 7.7 Jy) at the centre of 
the 3C196 field, rms of the residual is already low (1% of the rms 
of the original image) after calibrating with only 3C196 which con¬ 
tains 10% of the total flux of the field. Errors go down significantly 
when we include 15% of the total flux by adding another 4 sources 
in the model, but after that there is no rapid improvement. 


4 SIMULATION OF GALACTIC FOREGROUND 

So far we have considered leakages from unpolarized point sources 
into Stokes Q,U,V only, but, as mentioned before, our interest lies 
in the opposite case, i.e. leakage from polarization to total inten¬ 
sity. Compact radio sources are very weakly polarized and most 
of the point sources seen in polarization maps can be attributed to 
instrumental polarization and leakage. As at frequencies of tens to 
hundreds of MHz the polarized sky is dominated by Galactic dif¬ 
fuse synchrotron emission, we take real data of the 3C196 field 
observed by LOFAR, and create the simulated dataset using it as a 
sky model following the pipeline described in the next section. In 
these simulations, except for the one represented by Fig. |12| i, our 
Stokes I data contain only the noise leaked from Stokes Q, U. We 
do not add realistic noise to Stokes 1 until the final test because 
that would make the quantification of the intrinsic instrumental po¬ 
larization over the complete fe-space difficult, as the expected level 
of leakage is lower than the system noise. However, as a final test 
we add system noise to check the efficiency of a leakage removal 
technique. 


/-s-Pleakage removed by AW-projection [%] 
40 50 60 70 80 


Ci _ 1 



A8 (deg) 

Figure 9. Linear polarization leakage removed by AWImager as a percent¬ 
age of the leakage; both size and colour of the bubbles represent the same 
quantity. All sources in this simulated dataset had a Stokes I flux of 10 
Jy and their linear leakages were around O.I mjy. We see that up to 80% 
leakage could be removed using AWImager in this case. 


4.1 Simulation setup 

The general pipeline of the simulation of Galactic foreground is al¬ 
most same as that of the extragalactic foreground (boxes with solid 
borders in Fig. |^, but there are two major differences: here we 
simulate datasets for 161 spectral bands instead of just one, and ex¬ 
amine the leakages from Stokes Q,U to /, L rather than that from 
I to Q,U,V. The former enables us to examine the frequency be¬ 
haviour of the leakages through rotation measure synthesis and 3D 
power spectrum analysis, and the latter provides us with a realistic 
estimate of the amount of leakage into Stokes I to be expected in 
the current LOFAR-EoR observations of the 3C196 field. 

To make the 3C196 polarization sky model, we took Stokes 
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Figure 10. Faraday dispersion images of Galactic diffuse polarized emission within the central 4 degrees of the 3C196 field (top) and their leakages to Stokes 
I (middle) and V (bottom) caused by the LOFAR differential beams at the Faraday depths of 0 (left), +1 (middle) and +2 (right) rad/m^. The diameters of the 
inner and the outer circles are 2° and 3.8° (FWHM of LOFAR station beam at 150 MHz) respectively. The images have 480 X 480 pixels of 0.5' with a PSF 
of 3'. 


Q,U images for 161 subbands spanning 32 MHz centred at 150 
MHz that were produced from a single-night (8 hr) LOFAR obser¬ 
vation using the standard LOFAR calibration and imaging pipeline 
(e.g. see |Yatawatta et al.|2013[|Jelic et al.|2014| ). During the reduc¬ 
tion process, Dl-errors were removed and the data was also cor¬ 
rected for the element beam at the phase centre, thereby removing 
most of the instrumentally polarized point sources. We removed the 
remaining point sources by just masking them with noise so that 
diffuse emission dominates the image. The most significant sys¬ 
tematic errors that still remains in these images are the ionospheric 
Faraday rotation and the differential beam. A dataset with iono¬ 
spheric correction implemented is not necessary for our case as we 
are not concerned with analysis of the real data here, but only with 
the fraction of leakage; thus, any reasonable input model would 
serve our purpose. Also note that we are applying the ‘model’ dif¬ 


ferential beam to an image that already has the ‘true’ differential 
beam in it. This cannot be avoided as direction dependent calibra¬ 
tion or differential beam correction are yet to be done in this ob¬ 
serving window, but DD-correction would not bring any dramatic 
change in the final results that we want to produce, as the polariza¬ 
tion maps are dominated by diffuse emission and differential beam 
is only an 1 % effect. 

Stokes I and V in the model images were put to zero so that 
after applying the beam, they contain only leakages from Q, U. The 
following steps were performed to produce the final results. 

(i) DDE-corrupted visibilities at different frequencies are simu¬ 
lated using AWImager, as a prediction using BBS would currently 
take too much time. Here AWImager, in effect, carries out the for¬ 
ward transform of the major cycle and stops. 

(ii) Images from the simulated visibilities are produced using 
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Table 3. Setup of the simulation of Galactic foreground. 


Baselines used for simulated observation 

up to 3 kX 

Number of spectral subbands 

161 

Number of channels in each subband 

1 

Width of the channels / frequency resolution, 5u 

0.19 MHz 

Central frequency of the observing band 

150 MHz 

Total bandwidth, Aiy 

32 MHz 

Total observation time 

8 hr 

Integration time / time resolution 

10 s 

Baseline cut for imaging and PS estimation 

30 - 800 A 

Angular resolution (PSF) of the images 

4.3 arcmin 

Number of pixels in the images 

480 X 480 

Size of each pixel 

0.5 arcmin 

Maximum detectable Faraday depth, y/S/SX^ 

160rad/m2 

Largest resolvable structure in Faraday depth, 

0.96 rad/m^ 

Resolution in Faraday depth space, 2\/3/AA^ 

1 rad/m^ 

Minimum and maximum [Mpc“^] 

0.02-0.53 

Minimum and maximum fc II [Mpc~^] 

0.011-1.85 


CASA. Different parameters of the input model images and the final 
CASA images were kept the same; for details see table[^ 

(iii) We make 4 image-cubes by combining the images for 4 
Stokes parameters and also convert the fluxes in Jy to intensities 
in temperature following Eq.|27| 

(iv) To analyse Faraday structure of the leakage of polarized 
emission, RM-synthesis is performed on the cubes according to the 
formalism of section [23] re suiting in 3 ‘dirty’ (without deconvolv¬ 
ing the RMSF) RM-cubes: T’(<f>), Fi{^) and T'y (4?)[^ 

(v) Cylindrical and spherical 3D power spectra for all Stokes 
parameters are calculated from the image-cubes according to the 
formalism described in section l?. 6 . 2 l 


4.2 Results 

We first show the results of RM-synthesis of both polarization and 
leakage, but with a focus on the leakages, and then present the 3D 
power spectra produced from the image-cubes. 


4.2.1 RM synthesis 

The polarization RM-cube, F{^) basically represents the real data 
that we took as our input model, as leakage from Q to U and vice 
versa will be very small compared to their brightness. The maxi¬ 
mum intensity in the cube is ~ 5 K which is seen at $ = -1-1.2 
rad/m^ and the brightest structures are located within the Faraday 
depths of -1.5 and -1-5.0. Three slices of the F(^) cube at T>=0,1,2 
rad/m^ are shown in the top 3 panels of Fig. and the diffuse 
Galactic emission is prominent in all of them, but increases toward 
<E> = -1-1 rad/m^. For a detailed analysis of polarized emission in 
the 3C196 window seen by LOFAR, we refer the reader to Jelic 
et al. (in preparation). The corresponding slices of Fi{^) leakage 
are shown in the middle panels of the figure and here we see that 
the highest leakages appear at 4? = 0 and their peak is ~ 10 mK. 
As differential beams vary slowly with frequency (e.g. see Fig.[^, 
the leakages caused by them are a smooth function of frequency, 
thereby making them localized around 4? = 0 in RM space. This 


The RM-synthesis code of Michiel Brentjens was used for this purpose. 



Faraday depth, [rad m 

Figure 11. Four lines of sight along Faraday depth (Faraday spectmm) for 
4 bright pixels in the Faraday dispersion images T’(4>) (top), leakage into 
Fj)#) (middle) and leakage into Fy{^) (bottom). The pixels were chosen 
according to their intensity in F/(4>) and their RA and DEC are shown in 
the legend. 


property can be utilized to correct the effects of leakage, but per¬ 
forming a realistic leakage removal is beyond the scope of this pa¬ 
per (e.g. see |Geil et al.|201 1) . However, in section [43| we will show 
the results of a correction that does not take the differential beam 
into account. Another aspect of these images can be seen by focus¬ 
ing on the central 2 degrees; for example, although F(4? = 0) is 
highest in the central part, the corresponding T/(4> = 0 ) leakage 
is still much lower than that of the outer region (between the inner 
and the outer circle). This is expected as leakage terms of the beam 
Mueller matrix increase toward the outskirts. 

The 3 bottom panels of Fig.[T^show the same Faraday disper¬ 
sion images for leakages into Stokes V and they are much lower 
than the corresponding leakages into Stokes /—so much lower 
that they are dominated by leaked noise. This can be understood in 
terms of the differential beam of Fig. 0 - the M 42 and M 43 com¬ 
ponents of this matrix are responsible for leaking Q,U to H and 
we see that they are 2-3 orders of magnitude lower than the com¬ 
ponents responsible for leakage into Stokes I, i.e. M 12 and M 13 . 

Behaviour of the leakage in Faraday space can be seen more 
clearly in Fig. where we show 4 lines of sight along Faraday 
depth (Faraday spectrum) for 4 bright pixels in F(4?) (top), T/(4?) 
(middle) and TV(4?) (bottom). The bright pixels were chosen in 
f’/(4>) and then the corresponding pixels were found in ^(4?) and 
TV (4?). The fact that instrumental polarization and leakage appears 
at 4? = 0 in a Faraday spectrum, convolved with the RMSF, is 
evident from the middle panel of the figure. It is not so evident in 
TV (4?) due to the dominance of leaked noise. 
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Figure 12. Top'. Cylindrically averaged 3D power spectra (PS) of the polarized emission (V) within the central 4° of the 3C196 field (a) and its leakages into 
Stokes V (b) and I (c). Panel d shows the ratio between panel c and the coiresponding PS of the 21-cm differential brightness temperature for the fiducial 
model of |Mesinger et al.H2011^ at z = 9. The contours are drawn where this ratio is 1 for both the normal leakage and the leakage reduced by 70%. Bottom'. 
(e) Square-root of the ratio between the panels c and a expressed as a percentage. The g and / panels represent the same quantity as that of the c and e panels 
respectively, but for the case when 4 foreground components were removed from the leakage into Stokes I. Panel h shows the same percentage as the e and / 
panels, but when 60 mK noise was added to the /-leakage before foreground removal. The subscript R stands for GMCA residual. 


4.2.2 3D Power spectra 

A 3D power spectrum analysis of the DDE-corrupted image-cubes 
would be most interesting, as this would allow us to calculate the 
amount of polarized Galactic foreground leaked into a possible 
‘EoR window’ (a region in 3D Eourier space where the EoR signal 
is taken to be least contaminated) of LOFAR; for an example of an 
EoR window, see Fig. 1 of |Dillon et al.| ( (20T4l l that was made using 
the instrumental parameters of MWA. In the top panels of Fig.|12| 
we show the 3D cylindrical power spectra of the beam-corrupted 
polarized emission ("P), its leakages into Stokes / and V, and the 
ratio between the power spectra of the /-leakage and that of the 21- 
cm differential brightness temperature STi, of the fiducial model of 
|Mesinger et al.| ( [2011| >. The plots show the power that lies within a 
given fej_, fc|| bin in units of [mK]^ or as a ratio. 

The P spectrum (panel a) exhibits the same characteristics 
that one would expect based on the behaviour of the polarized 
emission in RM-space (described in the previous section). As in 
RM-space the brightest polarized emission were found near $ = 1 
rad/m^, so here the power is high at low A:|| (< 0.1 Mpc“^). Some 
additional power is seen in a wedge-shaped region at fey >0.1 and 
fcx >0.1 which can be attributed to the frequential unsmoothing 
of the intrinsically smooth polarized foreground by the frequency- 


varying PSF, and the extra power at high , k± is due to noise. 
At fcx < 0.04 power is very low, as expected, and it reaches its 
maximum at around 0.3 Mpc“'^. The maximum power is around 
4.4 X 10® [mK]^ which is found at the highest , fc||. 

The /-leakage spectrum (c) looks very similar to the P spec¬ 
trum, and the leakage power reaches up to ~ 5.5 [mK]^. At 
k± < 0.1 and fe|| >0.1 leakage power is 2-3 orders of magni¬ 
tude lower than the maximum. In order to see if / is just a scaled 
down version of P, we calculate the ratio s/I/P as a percentage of 
P, shown in panel e, which gives an estimate of the percentage of 
rms leakage at different k± and fey. Evidently, at fcy < 0.06 leak¬ 
age rms is 0.2%-0.3% of the polarization rms, and can go as high 
as 0.4% at high fc|| where noise leaked from P to 1 dominates. 

The V -leakage spectrum is shown in panel 6, and its level is 
much lower, the peak being around 0.06 [mK]^. The region at high 
fcy and high k± is dominated by noise, as signal-to-noise ratio is 
lower for longer baselines. By comparing this spectrum with that 
of P, it can be seen that the rms of the V -leakage at low fey is only 
~ 0.003% of the rms of the polarized emission which means that 
the uncorrected V leakage is negligible compared to the current 
noise levels in the EoR experiments within a FoV of 4°. 

As noted in the previous section, leakage is lower near the 
centre of the field (see Fig.|10|l. In order to quantify the associated 
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Figure 13. Spherically averaged 3D power spectra of the polarized emission 
within the central 4° of the 3C196 field (top solid), its leakages into Stokes 
I (middle solid) and V (bottom solid) caused by the LOFAR model beam, 
and the PS of the 21-cm differential brightness temperature at 2 = 9 
(solid with circles) for the fiducial model of|Mesinger et al.H2011). 


decrease in power, we calculated the power spectra of I, P and 
\P~JP within the inner 3° of the field. We found that, in this case, 
maximum leakage into I at low fc|| is ~ 4.9 [mK]^ and \fTJP ~ 
0 .2% at fc|| < 0.06 which is lower than the level of leakage within 
the inner 4°. 

To see how P —>■ 7 leakage affects the EoR signal, we took 
the 3D spherical power spectrum of the fiducial model of 21-cm 
differential brightness temperature at 2 = 9 from |Mesinger| 
|et al.| ( [2bll| > and calculated the corresponding cylindrical power 
spectrum as Aii(fci,fc||) = k\k\\A%i{k)/\2,{k]_ -|-The 
spherical PS is plotted in Fig.[T^along with the spectra of P, I and 
V, and the figure clearly shows that the EoR signal power is higher 
than the /-leakage at fe < 0.3 Mpc~^, and can be 2 orders of mag¬ 
nitude higher at the lowest scales. The ratio between fc|| 

and A|(A:^, fc|l) is plotted in Fig. ! 


12 


1 where the contours are drawn 


at Af/A^i = 1 for the normal case and the case when 70% of the 
leakage had been removed. Evidently, there is an ‘EoR window’ 
above the PSF-induced wedge and below fcn ~ 0.5 Mpc and 
the window extends up to fey ~ 1 Mpc“^ when 70% leakage is 
removed. 


4.3 Polarization leakage removal 

As |Jelic et al.| ( |20f0l section 7.2) have discussed at length, a leak¬ 
age of the polarized foreground into total intensity will be a major 
obstacle in detecting the EoR signal if (1) the level of leakage is 
comparable to the intensity of the EoR signal, and/or (2) frequency 
spectrum of the leakage mimics that of the signal. Fortunately, in 
the 3C196 field the latter is not the case, as we have seen that there 
is no significant polarization at high Faraday depths, or, equiva¬ 
lently, at high /c||. Flowever, the power of leakage could be compa¬ 
rable to that of the signal at high A:|| and hence leakage needs to be 
removed with sufficient accuracy to extend the EoR window. 

There are many methods for removing foregrounds from 
Stokes /; some assume spectral smoothness of the foreground and 
try to fit it out using polynomials, while others do not assume any¬ 
thing and hence are called ‘blind’ or non-parametric methods (for 
a list see |Chapman et al.|2014^ . The best way to remove the leak¬ 
age contribution of the foreground is, of course, to use the time- 


frequency-baseline dependent Mueller matrices during calibration 
and/or imaging to produce beam- and leakage-corrected images. 
Another potential way is to correct them in the Faraday dispersion 
images, i.e. correcting the F/(T>) using information from F{^) as 
demonstrated by |Geil et al.H20lT| see section 6.2 and Fig. 6). For 
leakages as smooth as in the field of 3C196, simply filtering the 
//($) for <!> ~ 0 could be another potential solution. Flowever, 
testing these methods is beyond the scope of this paper, and here we 
use a non-parametric foreground removal method, called GMCA 
(generalized morphological component analysis;|Bobin et al.|2007| 
|2008a|b^ , that has been shown to be able to remove foregrounds 
from simulated LOFAR-EoR data with high accuracy jChapman et| 
|al.|2013t . 

If a signal is represented as X = AS + N where S is the 
foreground to be extracted, N is noise and A is the mixing matrix, 
then GMCA tries to calculate a mixing matrix for which S is spars¬ 
est (have the least number of non-zero wavelet coefficients) in the 
wavelet domain. For details of the algorithm we refer the readers 
to |Chapman et al.| ( |2013| l. We run GMCA on the Stokes /-leakage 
cube to extract and subtract 4 components of the leaked foreground, 
as this number has been shown to yield good results ^Chapman ] 
|et al.|20i4} , and produce 3D cylindrical power spectrum from the 
residual cube which is shown in Fig.|12^. It clearly shows that the 
power of the smooth foregrounds at low k\\ (< 0.06) has been re¬ 
duced by almost two orders of magnitude by GMCA; compare it 
with the input /-leakage spectra of panel c that is plotted on the 
same scale. On the other hand, everything above fc|| =0.1 has been 
kept completely untouched due to low SNR—where S is the fore¬ 
ground and N is the noise including the cosmic signal—as GMCA 
cannot produce reliable model for the foregrounds when the SNR 
is low. In panel /, we plot the ratio of the GMCA residual PS and 
the polarization PS which shows that after GMCA subtraction, rms 
residual leakage at A:|| < 0.1 is around 0.1% of the polarized in¬ 
tensity. However, the EoR signal could also be removed along with 
the foreground in this case as there was no noise in Stokes I except 
for a very low level of noise leaked from Q, U. 

To see how additive noise affects the removal of leakage, we 
add 60 mK (rms) noise, which should be reached after 600 hours 
of integration using LOFAR ( [Chapman et al.|2014| ), to the Stokes I 
leakage maps at all frequencies. The noise was added to the visibil¬ 
ities and a new image cube was produced from the noisy visibili¬ 
ties. We run GMCA on the noisy /-leakage cube and produce a 3D 
cylindrical PS from the residual and take the square-root of the ratio 
of this PS with respect to the /-leakage which is shown as a percent¬ 
age in Fig. We see that almost no leakage has been removed 
in this case, not even in the relatively high SNR region at low k±. 
Therefore, we conclude that in case of such levels of noise, either 
a different strategy should be taken to remove foreground-leakage, 
or the leakage dominated region (where the leakage is more than 
the EoR signal) should be avoided to some extent (see jChapmanj 
jet al.j [2014| ( for a discussion on the relative merits of foreground 
removal and avoidance.) 


5 SUMMARY AND CONCLUSIONS 

This paper presents a first step in the analysis of the systematic er¬ 
rors of a radio interferometer, with a focus on polarization leakage, 
by simulating the LOFAR observations of both compact and diffuse 
emission in the presence of direction independent and direction de¬ 
pendent errors which are treated separately. We have revisited the 
measurement equation of a radio interferometer and modelled the 
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direction independent (DI) and direction dependent (DD) errors as 
2x2 Jones matrices and the corresponding 4x4 Mueller matrices 
have been used to show the polarization properties of the instru¬ 
ment. The full polarization DD-Mueller matrix (Fig. describing 
the time-frequency-direction dependent behaviour (e.g. see Fig. 
and of the response of a baseline of LOFAR, created by two 
stations (Fig. [^, has been presented to be a DD equivalent of the 
DI-Mueller matrix of Eq.[T^ 

We have simulated an observation with Dl-errors by assum¬ 
ing them to be random at every timestep and the rms of the ran¬ 
dom numbers drawn from a Gaussian distribution with zero mean 
is dubbed the ‘rms Dl-error’. We find that self-calibration can solve 
for these errors to an extremely high accuracy if the sky model is 
perfect as, in that case, the information provided by an interfer¬ 
ometer will be highly redundant. For an rms Dl-error of 10“^, the 
selfcal error is less than 0.002% and the corresponding error in the 
rms of the resulting residual image is less than 0.005% (Fig.[^. 

The only DD-error that we have simulated is the differen¬ 
tial (normalized with respect to the phase centre) station beam of 
LOFAR. We simulated the LOFAR observations of extragalactic 
unpolarized point sources in the 3C196 observing window of the 
LOFAR-EoR experiment including the DD-errors and estimated 
the flux and position errors due to self-calibration with incom¬ 
plete sky models and the percentage of / —>■ (Q, U) leakage of 
the brightest sources (see Fig.|^. We see that the errors go down 
significantly as the sky model is improved. Flowever, calibrating 
with only unpolarized sources has its limitations, e.g. the unitary 
ambiguity ( [Wijnholds, et al.|2012[|Carozzi|2014) . There is no plan 
for using polarized sources in calibrating LOFAR EoR data until 
now, as there are very few intrinsically polarized point sources in 
the data, and the polarized emission is dominated by diffuse emis¬ 
sion (e.g. see |Yatawatta et al.|[20T3l |Jelic et ar]|2014^ . We test a 
possible strategy of correcting the DD errors from point sources 
using AWImager with an unrealistic, exaggerated sky model and 
see that AWImager can remove up to 80% of the leakage from 
Stokes Q, U, but a more elaborate testing of this algorithm with 
realistic sky models has to be done to reach any final conclusion. 


To predict the level of polarization leakage in the Stokes I 
images of the 3C196 field, we took the real LOFAR observations 
of Galactic diffuse polarized emission in this field and created an 
unreal sky model where I = V = 0 to quantify the leakages 
from Q, U lo 1,V caused by the DD errors. An RM-synthesis of 
the DDE-corrupted V image cubes showed that in this particular 
field polarization peaks within the Faraday depths (4>) of -1 and 
-(-5 rad/m^. From the effective Stokes I Faraday dispersion images 
we saw that polarization leakage is localized around = 0 (Fig. 
[TTJ, as DD-errors do not have any rapid variation along frequency. 
Maximum leakage was found to be around 15 mK which could be 
comparable to the EoR signal (Fig.|lQ|>. 


To understand the level of leakage contaminant in the ‘EoR 
window’ of the instrumental fc-space, we calculated the cylindri- 
cally and the spherically averaged 3D power spectra (PS) of I, P, V 
cubes. The P spectrum shows characteristic smooth polarized fore¬ 
grounds at low fc|| (Fig. 12 I and the /-leakage spectrum looks very 
similar to this. From the power ratio, \/I/P we showed that the 
percentage of rms leakage over the k±, space varies by a factor 
of 2 and ranges from 0.2% to 0.4%. We compared the /-leakage 
with the 3D PS of the expected 21-cm differential brightness tem¬ 
perature at 2 : = 9 simulated by |Mesinger et alT| ( |2011^ and saw 
that the region above the PSF-induced wedge and below fc|| ~ 0.5 
Mpc“^ is dominated by the cosmic signal (Fig. 


12 


i) and hence 


defines a potential ‘EoR window’, and the window expands sub¬ 
stantially after removing 70% of the leakage. 

As the /-leakage do not mimic the EoR signal in this case, 
we tried to remove it using GMCA which is being used to remove 
diffuse foreground from the LOFAR-EoR data. Erom the 3D PS of 
the residual left after the removal of foreground leakage compo¬ 
nents by GMCA, we saw that (Pig. 12 ',g) at k^ < 0.1, i.e. in the 


high SNR regime, GMCA could reduce the leakage by up to two 
orders of magnitude while the region above that scale was left com¬ 
pletely untouched. For a more realistic analysis, we added 60 mK 
noise to the Stokes I leakage maps, reran GMCA on it and saw that 
(Fig .|12}i) in this case almost no leakage was removed, not even in 
the relatively high SNR region. 

Antennas for the future arrays like SKA, that have EoR detec¬ 
tion as one of the main scientific objectives, are being designed in 
such a way that their polarimetric performance is good enough to be 
able to minimize the effects of polarization leakage (de Lera Acedo; 
private communication). A recently proposed figure of merit for 
quantifying the polarimetric performance is the intrinsic cross¬ 
polarization ratio (IXR) which, in Mueller formalism, can be di¬ 
rectly related to the instrumental polarization jCarozzi & Woan] 
|201 1[ eq. 23). Our LOFAR results show an instrumental polariza¬ 
tion of around 0.3% (Fig.|12^; ignoring V ^ I leakage) within the 
FWHM of the nominal station beams, i.e. within a FoV of ~ 4°. 
This corresponds to an IXRm {Mueller IXR) of 25 dB, or equiva¬ 
lently an IXRj {Jones IXR, see eq. 25 of |Carozzi & Woan|2011| ) 
of 56 dB, and if the leakage can be reduced by 70%, IXRm will 
improve to 35 dB. Therefore, we can say that if SKA has a min¬ 
imum IXRm of 25 dB within the central ~ 4° of its nominal sta¬ 
tion beams, then even a modest polarimetric calibration (~ 70% 
leakage removal) will ensure that the polarization leakage remains 
well below the expected EoR signal at the scales of 0.02-1 Mpc”'^. 
However, if the IXRm is lower within a FoV of 4 °, more leakage 
needs to be removed to reach the same level as before in relation 
to the EoR signal in the power spectra, e.g. if the IXRm is 20 dB, 
91% leakage has to be removed, and if it is 15 dB, 97% has to be 
removed. 

The major conclusions of this paper are the following. 


(i) Two properties of the polarization leakage can be utilized for 
its removal in this specific case: it appears around a Earaday depth 
of 0 rad/m^ in RM-space and the overall variation of the rms of the 
fractional leakage in the instrumental fe-space is less than a factor 
of 2. 

(ii) In the cylindrically averaged 3D power spectra, a clear ‘EoR 
window’ can be defined in terms of polarization leakage above the 
wedge and below A:|| ~ 0.5 Mpc“^. Within this window, the EoR 
signal dominates the polarization leakage and the window takes up 
the whole fc-space at k\\ < 1 after removing 70% of the leakage. 

(iii) A DDE-blind foreground removal method like GMCA is 
not ideal for removing leakage of diffuse polarized emission, as the 
level of leakage is lower than the current noise level in the LOFAR 
observations. 
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